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Low-dimensionality and predictability of solar wind and global 
magnetosphere during magnetic storms 

T. Zivkovic, 1 and K. Rypdal, 1 

Abstract. The storm index SYM-H, the solar wind velocity v, and interplanetary mag- 
netic field B z show no signatures of low-dimensional dynamics in quiet periods, but tests 
for determinism in the time series indicate that SYM-H exhibits a significant low-dimensional 
component during storm time, suggesting that self-organization takes place during mag- 
netic storms. Even though our analysis yields no discernible change in determinism dur- 
ing magnetic storms for the solar wind parameters, there are significant enhancement 
of the predictability and exponents measuring persistence. Thus, magnetic storms are 
typically preceded by an increase in the persistence of the solar wind dynamics, and this 
increase is also present in the magnetospheric response to the solar wind. 



1. Introduction 

Under the influence of the solar wind, the magnetosphere 
resides in a complex, non-equilibrium state. The plasma 
particles have non-Maxwellian velocity distribution, MHD 
turbulence is present everywhere, and intermittent energy 
transport known as bursty-bulk flows occurs as well [An- 
gelopoulos et al., 1999]. The magnetospheric response to 
particular solar events constitutes an essential aspect of 
space weather while the response to solar variability in gen- 
eral is often referred to as space climate [Watkins, 2002]. 
Theoretical approaches to space climate involve concepts 
and methods from stochastic processes, nonlinear dynamics 
and chaos, turbulence, self-organized criticality, and phase 
transitions. 

Self-organization can lead to low-dimensional behavior in 
the magnetosphere [Klimas et al, 1996; Vassiliadis et al, 
1990; Sharma et al., 1993]. However, power-law dependence 
observed in the Fourier spectra of the auroral electrojet 
(AE) index is a typical signature of high dimensional col- 
ored noise indicating multi-scale dynamics of the magneto- 
sphere. In order to reconcile low-dimensional, deterministic 
behavior with high-dimensionality, Chang [1998] proposed 
that a high-dimensional system near self-organized critical- 
ity (SOC) [Bak et al, 1987] can be characterized by a few 
parameters whose evolution is governed by a small number 
of nonlinear equations. Some magnetospheric models, like 
the one presented in Chapman et al. [1998], are based on the 
SOC-concept. Here a system tunes itself to criticality and 
the energy transport across scales is mediated by avalanches 
which are power-law distributed in size and duration. 

On the other hand, it was suggested in Sitnov et al. 
[2001] that the substorm dynamics can be described as a 
non-equilibrium phase transition; i.e. as a system tuned ex- 
ternally to criticality. Here, a power-law relation is given, 
with characteristic exponent close to the input-output crit- 
ical exponent in a second-order phase transition. In fact, it 
is claimed in Sharma et al. [2003] that the global features of 
the magnetosphere correspond to a first order phase transi- 
tion whereas multi-scale processes correspond to the second- 
order phase transitions. 

The existence of metastable states in the magnetosphere, 
where intermittent signatures might be due to dynamical 
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phase transitions among these states, was suggested by Con- 
solini and Chang [2001], and forced and/or self-organized 
criticality (FSOC) induced by the solar wind was introduced 
as a conceptual description of magnetospheric dynamics. 
The concept of intermittent criticality was suggested by Bal- 
asis et al. [2006] who asserted that during intense magnetic 
storms the system develops long-range correlations, which 
further indicates a transition from a less orderly to a more 
orderly state. Here, substorms might be the agents by which 
longer correlations are established. This concept implies a 
time-dependent variation in the activity as the critical point 
is approached, in contrast to SOC. 

In the present paper we investigate determinism and pre- 
dictability of observables characterizing the state of the 
magnetosphere during geomagnetic storms as well as dur- 
ing its quiet condition, but the emphasis is on the evolu- 
tion of these properties over the course of major magnetic 
storms. The measure of determinism employed here in- 
creases if the system dynamics is dominated by modes gov- 
erned by low-dimensional dynamics. Hence, the determin- 
ism in most cases is a measure of low-dimensionality. For a 
low-dimensional, chaotic system the predictability measure 
increases when the largest Lyapunov exponent increases, 
and hence it is really a measure of un-predictability. For 
a high-dimensional or stochastic system it is related to the 
degree of persistence in time series representing the dynam- 
ics. High persistence means high predictability. 

One of the most useful data tools for probing the mag- 
netosphere during substorm conditions is the AE minute 
index which is defined as the difference between the AU in- 
dex, which measures the eastward electrojet current in the 
auroral zone, and the AL index, which measures the west- 
ward electrojet current, and is usually derived from 12 mag- 
netometers positioned under the auroral oval [Davies and 
Sugiura, 1966]. The auroral electrojet, however, does not 
respond strongly to the specific modifications of the mag- 
netosphere that occur during magnetic storms. A typical 
storm characteristic, however, is a change in the intensity of 
the symmetric part of the ring current that encircles Earth 
at altitudes ranging from about 3 to 8 Earth radii, and is 
proportional to the total energy in the drifting particles that 
form this current system [Gonzalez et al., 1994]. The indices 
D st and SYM-H indices are both designed for the study of 
storm dynamics. These indices contain contribution from 
the magnetopause current, the partial and symmetric ring 
current, the substorm current wedge, the magnetotail cur- 
rents, and induced currents on the Earth's surface. They are 
derived from similar data sources, but SYM-H has the dis- 
tinct advantage of having 1-min time resolution compared 
to the 1-hour time resolution of D st . Wanliss [2006] has 
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recommended that the SYM-H index be used as a de facto 
high-resolution D st index. The analysis of these indices are 
central to this study. We particularly focus on SYM-H and 
SYM-H* which is derived from the SYM-H when the con- 
tribution of the magnetopause current is excluded. 

The typical magnetic storm consists of the initial phase, 
when the horizontal magnetic field suddenly increases and 
stays elevated for several hours, the main phase where this 
component is depressed for one to several hours, and the re- 
covery phase which also lasts several hours. The initial phase 
has been associated with northward directed IMF (little en- 
ergy enters the magnetosphere) , but it has been discovered 
that this phase is not essential for the storm to occur Akasofu 
[1965]. In order to define a storm, we follow the approach of 
Loewe and Prolss [1997], where the D st minimum is a com- 
mon reference epoch, the main-phase decrease is sufficiently 
steep, and the recovery phase is also defined. 

1.1. Data acquisition 

The SYM-H index data are downloaded from World Data 
Center, with 1-min resolution. We also use minute data 
for the interplanetary magnetic field (IMF) component B z , 
minute data for the solar wind bulk velocity v along the 
Sun-Earth axis, as well as flow pressure which is given in nT. 
These data are retrieved from the OMNI satellite database 
and are given in the GSE coordinate system. Gaps of miss- 
ing data in B z , v and flow pressure are linearly interpolated 
from the data which are not missing, while SYM-H data are 
analyzed for the entire period. The same result for the B z 
and v is obtained when gaps of missing data are excluded 
from the analysis. 

Data for the period from January 2000 till December 2005 
is used to compute general properties of the magnetosphere. 
In order to analyze storm conditions all the indices are ana- 
lyzed during ten intense magnetic storms. Analyzed storms 
occurred on 6 April 2000, 15 July 2000, 12 August 2000, 31 
March 2001, 21 October 2001, 28 October 2001, 6 November 
2001, 7 September 2002, 29 October 2003, and 20 Novem- 
ber 2003. These storms are characterized with D st minimum 
which is in the range between -150 nT to -422 nT. 

The remainder of the paper is organized as follows: sec- 
tion 2 describes the data analysis methods employed. Sec- 
tion 3 presents analysis results discerning general statistical 
scaling properties of global magnetospheric dynamics using 
minute data over several years and data generated by a nu- 
merical model which produces realizations of a fractional 
Ornstein-Uhlenbeck (fO-U) process. In particular we study 
how determinism and predictability of the geomagnetic and 
solar wind observables change over the course of magnetic 
storms. Section 4 is reserved for discussion of results and 
section 5 for conclusions. 

2. Methods 

2.1. Recurrence-plot analysis 

The recurrence plot is a powerful tool for the visualiza- 
tion of recurrences of phase-space trajectories. It is very 
useful since it can be applied to non-stationary as well as 
short time series [Eckmann el al, 1987], and this is the na- 
ture of data we use to explore magnetic storms. Prior to 
constructing a recurrence plot the common procedure is to 
reconstruct phase space from the time-series x(t) of length 
N by time-delay embedding [Takens, 1981]. 

Suppose the physical system at hand is a deterministic 
dynamical system describing the evolution of a state vector 
z(t) in a phase space of dimension p, i.e. z evolves according 
to an autonomous system of 1st order ordinary differential 
equations; 
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and that an observed time series x(t) is generated by the 
measurement function g : TV — > 1Z, 

x(t) = g(z(t)). (2) 

Assume that the dynamics takes place on an invariant set 
(an attractor) A C TZ P in phase space, and that this set has 
box-counting fractal dimension d. Since the dynamical sys- 
tem uniquely defines the entire phase-space trajectory once 
the state z(t) at a particular time t is given, we can define 
uniquely an m-dimensional measurement function, 

g : A K m , g(z) = (x{t), x(t + T),...,x(t+(m- l)r)). 

(3) 

where the vector components are given by equation (2), and 
t is a time delay of our choice. If the invariant set A is 
compact (closed and bounded), g is a smooth function and 
m > 2d, the map given by equation (3) is a topological em- 
bedding (a one-to-one continuous map) between A and lZ m . 
The condition m > 2d can be thought of as a condition for 
the image g(.4) not to intersect itself, i.e. to avoid that two 
different states on the attractor A are mapped to the same 
point in the m-dimensional embedding space lZ m . If such 
an embedding is achieved, the trajectory x(t) = g(z) (where 
g(z) is given by equation (3)) in the embedding space is a 
complete mathematical representation of the dynamics on 
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Figure 1. B z during quiet condition on September 5, 
2001. a) B z time series, b) Recurrence plot of the time 
series shown in (a). 
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the attractor. Note that the dimension p of the original 
phase space is irrelevant for the reconstruction of the em- 
bedding space. The important thing is the dimension d of 
the invariant set A on which the dynamics unfolds. 

There are practical constraints on useful choices of the 
time delay r. If r is much smaller than the autocorrelation 
time the image of A becomes essentially one-dimensional. If 
r is much larger than the autocorrelation time, noise may de- 
stroy the deterministic connection between the components 
of x(t), such that our assumption that z(t) determines x(t) 
will fail in practice. A common choice of r has been the 
first minimum of the autocorrelation function, but it has 
been shown that better results are achieved by selecting the 
time delay as the first minimum in the average mutual in- 
formation function, which can be percieved as a nonlinear 
autocorrelation function [Abarbanel, 1996]. Here we use the 
average mutual information function to calculate the value 
of r. 

The recurrence-plot analysis deals with the trajectories in 
the embedding space. If the original time series x(t) has N 
elements, we have a time series of N — (m — l)r vectors x(t) 
for t = 1,2, . . . , N — (m — l)r. This time series constitutes 
the trajectory in the reconstructed embedding space. 

The next step is to construct a [N—(m— l)r] x [(N— (m— 
1)t] matrix Rij consisting of elements and 1. The matrix 
element is 1 if the distance is ||x; — x 3 -|| < e in the 

reconstructed space, and otherwise it is 0. The recurrence 
plot is simply a plot where the points for which the 
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Figure 2. B z during the strong storm on 6 April 2000. 
a) B z time series, b) Recurrence plot of the time series 
shown in (a). 
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Figure 3. a) Intrinsic mode functions obtained by EMD 
for B z for the magnetic storm on 6 April 2000, b) D st for 
the same event. 



corresponding matrix element is 1 is marked by a dot. For a 
deterministic system the radius e is typically chosen as 10% 
of the diameter of the reconstructed attractor, but varies for 
different sets of data. For a non-stationary stochastic pro- 
cess like a Brownian motion there is no bounded attractor 
for the dynamics, and the diameter is limited by the length 
of the data record. The first example of recurrence plot is 
shown in Figure 1, obtained from the B z when no storm is 
present on 5 September 2001. In Figure 2, the recurrence 
plot is shown for the B z for the strong storm on 6 April 
2000. In both cases, embedding dimension is m = 1 and 
e ~ 0.4, which corresponds to 10% of the data range. 

2.2. Empirical mode decomposition 

The empirical mode decomposition (EMD) method, de- 
veloped in Huang et al. [1998] is very useful on non- 
stationary and nonlinear time series. EMD method can give 
a change of frequency in any moment of time (instantaneous 
frequency) and a change of amplitude in the system. How- 
ever, in order to properly define instantaneous frequency, a 
time series should have the same number of zero crossings 
and extrema (or they can differ at most by one), and a lo- 
cal mean should be close to zero. The original time series 
usually does not have these characteristics and should be 
decomposed into intrinsic mode functions (IF) for which in- 
stantaneous frequency can be defined. Decomposition can 
be obtained through the so-called sifting process. This is an 
adaptive process derived from the data and can be briefly 
described as follows: All local maxima and minima in the 
time series s(t) are found, and all local maxima and min- 
ima are fitted by cubic spline and these fits define the upper 
(lower) envelope of the time series. Then the mean of the 
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upper and lower envelope m{t) is defined, and the difference 
between the time series and this mean represents the first 
IF, h(t) — s(t) — m(t), if instantaneous frequency can be ob- 
tained, defined by some stopping criterion. If not, the proce- 
dure is repeated (now starting from h(t) instead of s(t)) until 
the first IF is produced. Higher IFs are obtained by sub- 
tracting the first IF from the time series s(t) and the entire 
previously mentioned procedure is repeated until a residual, 
usually a monotonic function, is left. We use a stopping cri- 
terion defined by Rilling et al. [2003], where r)(t) < 81 on 
1 — 7 fraction of the IF, and r/(t) < 62, on the remaining 
fraction of the IF. Here r\ = m(t)/a(t), a(t) is the IF ampli- 
tude, and 7 = 0.05, 9\ = 0.05, and 82 = 0.5. By the above 
definitions, IFs are complete in the sense that their summa- 
tion gives the original time series: s(t) — ^ h(t) + R(t) 
where M is the number of IFs and R is a residual. In Figure 
3a we show the IFs from EMD performed on the IMF B z 
during a magnetic storm on 6 April 2000 (whose time series 
is plotted in Figure 2a), while in Figure 3b the D s t index 
for the same storm is shown. 

In order to study stochastic behavior of a time series 
by means of EMD analysis, we refer to Wu and Huang 
[2004] who studied characteristics of white noise using the 
EMD method. They derived for white noise the relation- 
ship log E m — — log T m , where E m and T m represents em- 
pirical variance and mean period for the m'th IF. Here, 

E m = (1/A0 J2t=i /l ( t ) 2 ' where h(t) is the m'th IF and T m 
is the ratio of the m'th IF length to the number of its zero 
crossings. Franzke [2009] analyzed telecommunication in- 
dices and noticed a resemblance to autoregressive processes 
of the first order AR(1), which are stochastic and linear pro- 
cesses. For such processes log E m — £logT m . For fractional 
Gaussian noise processes (H < 1) and fractional Brown- 
ian motions (H > 1) we have the connection £ = 2H — 2, 
where H is the Hurst exponent, as shown by Flandrin and 
Gongalves [2004]. A useful feature of the EMD analysis is 
the possibility of extraction of trends in the time series [ Wu 
et al, 2007], because the slowest IF components should often 
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be interpreted as trends. This is an advantage compared to 
the standard variogram or rescaled-range techniques [Reran, 
1994], whose estimation of the scaling exponents is biased 
by the trend. 

2.3. A test for determinism 

In this paper we employ a simple test for determinism, 
developed by Kaplan and Glass [1992], where the following 
hypothesis is tested: When a system is deterministic, the 
orientation of the trajectory (its tangent) is a function of 
the position in the phase space. Further, this means that 
the tangent vectors of a trajectory which recurs to the same 
small "box" in phase space, will have the same directions 
since these are uniquely determined by the position in phase 
space. On the other hand, trajectories in a stochastic sys- 
tem have directions which do not depend uniquely on the 
position and are equally probable in any direction. This 
test works only for continuous flows, and is not applicable 
to maps since consecutive points on the orbit may be very 
separated in the phase space. For flows, the trajectory orien- 
tation is defined by a vector of a unit length, whose direction 
is given by the displacement between the point where tra- 
jectory enters the box j to the point where the trajectory 
exits the same box. The displacement in m-dimensional 
embedding space is given from the time-delay embedding 
reconstruction: 

Ax(t) = [x(t+b)-x(t),x(t + T + b)-x(t + T),..., 

x(t + (m - 1)t + b) - x(t + (m - l)r)], (4) 

where b is the time the trajectory spends inside a box. The 
orientation vector for the fcth pass through box j is the unit 
vector Ufcj = Axk,j(t)/\Axk,j(t)\. The estimated averaged 
displacement vector in the box is 
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Figure 4. a) Time series representing one component 
of the numerical solution of the Lorenz system, b) Av- 
erage displacement vectors Vj in each box visited by a 
2-dimensional projection of the m = 3 -dimensional em- 
bedding space reconstructed from the time series in (a). 



Figure 5. a) Time series representing the numerical 
solution of the equation for the f-OU process, b) Av- 
erage displacement vectors Vj in each box visited by a 
2-dimensional projection of the m = 8 -dimensional em- 
bedding space reconstructed from the time series in (a). 
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where rij is the number of passes of the trajectory through 
box j. If the dynamics is deterministic, the embedding di- 
mension is sufficiently high, and in the limit of vanishingly 
small box size, the trajectory directions should be aligned 



and ^ = 1^1 



1. In the case of finite box size, Vj will not 



depend very much on the number of passes rij, and Vj will 
converge to 1 as nj — > oo. In contrast, for the trajectory 
of a random process, where the direction of the next step is 
completely independent of the past, Vj will decrease with rij 
as Vj ~ rij . In our analysis we will choose the linear box 
dimension equal to the mean distance a phase-space point 
moves in one time step and set 6=1 time step in equation 
(4). 

In Figure 4b we show displacement vectors Vj averaged 
over the passes through the box j, for a three-dimensional 
embedding of the Lorenz attractor, whose time series is 
shown in Figure 4a; in Figure 5b the same is shown for a 
random process, in this case a fractional Ornstein-Uhlenbeck 
(fO-U) process. These model systems will be used through- 
out this paper as archetypes of low-dimensional and stochas- 
tic systems, respectively. The Lorenz system has the form 



dx/dt — a(y — x) 
dy/dt — —xz + cx — y 
dz/dt — xy — bz, 



(6) 



with standard coefficient values a = 10, b — 8/3, and c = 28, 
which give rise to a chaotic flow. The fO-U process is de- 
scribed by the stochastic equation: 



dSt = \{ti-St)+adW t , 



(7) 



where dWt is a fractional Gaussian noise with Hurst expo- 
nent H [Beran, 1994]. The drift (A and /x) and diffusion 
(a) parameters are fitted by the least square regression to 
the time series of the SYM-H storm index. This will be 
explained in more detail in section 3.2. 

The degree of determinism of the dynamics can be as- 
sessed by exploring the dependence of Vj on rij . In practice, 
this can be done by computation of the averaged displace- 
ment vector 



{Vj)u j= r, 



(8) 



where the average is done over all boxes with same number n 
of trajectory passes. As shown in Kaplan and Glass [1993], 
the average displacement of n passes in m-dimensional phase 
space for the Brownian motion is 



,F[(m+l)/2] 



^n y m> r(m/2) 



(9) 



where F is the gamma function. The deviation in (Vj) be- 
tween a given time series and the Brownian motion can be 
characterized by a single number given by the weighted av- 
erage over all boxes of the quantity, 



A(r) 



1 



Rrn 



(10) 



where we have explicitly highlighted that the averaged dis- 
placement (Vj)(r) of the trajectory in the reconstructed 
phase space depends on the time-lag r. For a completely 
deterministic signal we have A = 1, and for a completely 
random signal A = 0. 

All systems described by the laws of classical (non- 
quantum) physics are deterministic in the sense that they 
are described by equations that have unique solutions if the 
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Figure 6. a) L n : square symbols are derived from nu- 
merical solutions of the Lorenz system, and triangles from 
these solutions after randomization of phases of Fourier 
coefficients, b) A(r): diamonds from Lorenz system, and 
triangles after randomization of phases. 



initial state is completely specified. In this sense it seems 
meaningless to provide tests for determinism. The test de- 
scribed in this section is really a test of low dimensionality. 
The test is performed by means of a time-delay embedding, 
for embedding dimension m up to a maximum value M, 
where M is limited by practical constraints. High M re- 
quires longer time series in order to achieve adequate statis- 
tics. A test that fails to characterize the system as determin- 
istic for m < M in reality only tells us that the embedding 
dimension is too small, i.e. the number of degrees of free- 
dom d of the system exceeds M/2. Such systems will in the 
following be characterized as random, or stochastic. 

In Figure 6a, we plot L n versus n for a time series gener- 
ated as a numerical solution of the Lorenz system. Here we 
use m = 3, & = 1 , t = 14 and the box size is of the order of 
average distance a phase-space point moves during one time- 
step. In the same plot we also show the same characteristic 
for the surrogate time series generated by randomizing the 
phases of the Fourier coefficients of the original time series. 
This procedure does not change the power spectrum or auto- 
covariance, but destroys correlation between phases due to 
nonlinear dynamics. For low-dimensional, nonlinear systems 
such randomization will change L n , as is demonstrated for 
the Lorenz system in Figure 6a. We also calculate A ver- 
sus r for these time series and plot the results in Figure 6b. 
Again, A(r) for the original and surrogate time series are 
significantly different. 

For the numerically generated fO-U process, where m = 
8, b = 1 and r = 20, we observe in Figure 7 that L n and 
A(t) for the original and surrogate time series do not differ, 
demonstrating that these quantities are insensitive to ran- 
domization of phases of Fourier coefficients if the process is 
generated by a linear stochastic equation. 
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One should pay attention to the nature of the exper- 
imental data used in the test of determinism. For low- 
dimensional data contaminated by low-amplitude noise or 
Brownian motions, the analysis results will depend on the 
box size, but the problem is solved by choosing it sufficiently 
large. For a low-dimensional system represented by an at- 
tractor of dimension d the results may also depend on the 
choice of embedding dimension m. The estimated determin- 
ism L n tends to increase with increasing m until it stabilizes 




Figure 7. a) L n : square symbols are derived from nu- 
merical solutions of the fO-U stochastic equation, and tri- 
angles from these solutions after randomization of phases 
of Fourier coefficients, b) A(r): squares from fO-U equa- 
tion, and triangles after randomization of phases. 



at L n w 1 as m approaches 2d. For a random signal there 
is no such dependence on embedding dimension, as demon- 
strated by example in Figure 8. Here we plot the determin- 
ism Z/3 (Z/„ when n — 3) versus embedding dimension m for 
the Lorenz and fO-U time series. For comparison we also 
plot this for transformed SYM-H, tSYM-H, during magnetic 
storm times (the transformation and reasons for it are ex- 
plained in section 3). It converges to a value less than 1, and 
for embedding dimensions higher than for the Lorenz time 
series. This indicates that this geomagnetic index during 
magnetic storms exhibit both a random and a deterministic 
component, and that the dimensionality of this component 
is higher than for the Lorenz system. 

2.4. A test for predictability 

In this subsection we develop an analysis which is based 
on the diagonal line structures of the recurrence plot. In our 
study we use the average inverse diagonal line length: 



r=(z- 1 )=^r 1 P(Z)/^P(0, 



(ii) 



where P(l) is a histogram over diagonal lengths: 



R 



i-W-l 



)(l-R. 



i + k,j + k- 



For a low-dimensional, chaotic deterministic system (for 
which the embedding dimension is sufficient to unfold the 
attractor) V is an analog to the largest Lyapunov exponent, 
and is a measure of the degree of unpredictability. 

For stochastic systems, the recurrence plots do not have 
identifiable diagonal lines, but rather consists of a pattern 
of dark rectangles of varying size, as observed in Figure 1. 
For embedding dimension m = 1 such a dark rectangle cor- 
responds to time intervals 7i = (ti,ti + Ati) on the hori- 
zontal axis and I2 ~ (£2,^2 + on the vertical axis, for 
which the signal x(t) is inside the same e-interval whenever 
t is included in either Ii or Ii- In this case the length of 
unbroken diagonal lines I is a characteristic measure of the 
linear size of the corresponding rectangle, and the PDF P(l) 
a measure of the distribution of residence times of the tra- 
jectory inside e-intervals. For selfsimilar stochastic processes 
such as fractional Brownian motions P(l) can be computed 
analytically, and V computed as function of the selfsmilar- 
ity exponent h. Since the residence time I inside an e-box 
increases as the smoothness of the trajectory increases (in- 
creasing h), we should find that T(h) is a monotonically 




'2 4 6 8 10 _ 10 l , , , i_ 

m 0.5 1 1.5 2 



time (min) 



Figure 8. £3 as a function of embedding dimension m 
for solution of Lorenz equations (triangles), fO-U process 
(squares), and tSYM-H (circles). 



Figure 9. a) Increments for the SYM-H index, b) In- 
crements for the transformed signal tSYM-H. 
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decreasing function of h. In section 3.4 we compute T(h) 
numerically for a synthetically generated fO-U process and 
thus demonstrate this relationship between F and h. Hence 
both r and h can serve as measures of predictability, but V is 
more general, because it is not restricted to selfsimilar pro- 
cesses or processes with stationary increments, and applies 
to low-dimensional chaotic as well as stochastic systems. 

3. Results 

In Rypdal and Rypdal [2010] it is shown that the fluctu- 
ation amplitude (or more precisely; the one-timestep incre- 
ment) Ay(t) of the AE index is on the average proportional 
to the instantaneous value y(t) of the index. This gives 
rise to a special kind of intermittency associated with mul- 
tiplicative noises, and leads to a non-stationary time series 
of increments. However, the time series Ay(t)/y(t) is sta- 
tionary, implying that the stochastic process x(t) = logy(t) 
has stationary increments. Thus, a signal with stationary 
increments, which still can exhibit a multifractal intermit- 
tency, can be constructed by considering the logarithm of 
the AE index. Similar properties pertain to the SYM-H 
index, although in these cases we have to add a constant 
ci before taking the logarithm, i.e. x(t) = log (ci + C2y(t)) 
has stationary increments. Using the procedure described 
in Rypdal and Rypdal [2010] the estimated coefficients are 
ci = 0.7725 and C2 = 0.0397. In Figure 9a we show the in- 
crements for the original SYM-H data, while in Figure 9b we 
show the increments for the transformed signal x(t), which 
in the following will be denoted tSYM-H. 

3.1. Scaling of storm- and solar wind parameters 

In this section, we employ EMD and variogram analy- 
sis to tSYM-H, IMF B z and solar wind flow speed v. The 
EMD analysis is used to compute intrinsic mode functions 
(IF) for time intervals of 50000 minutes using data for the 
entire period from January 2000 till December 2005. The 
empirical variance estimates E versus mean period T for 
each IF component in tSYM-H, B z , and v are shown as log- 
log plots in Figure 10a. In section 2.2 we mentioned that 
Flandrin and Gongalves [2004] has demonstrated that for 
fractional Gaussian noise the slope £ is equal to ( — 2H — 2, 
where H is the Hurst exponent. This estimate for the slope 
seems valid for our data as well, as is shown in the figure 
from comparison with the variogram, even though the time 
series on scales up to 10 4 minutes are non-stationary pro- 
cesses having the character of fractional Brownian motions 
[Rypdal and Rypdal, 2010]. The results from the two dif- 
ferent methods shown in Figures 10a and 10b are roughly 
consistent, using the relations h = H — 1 and £ = 2H — 2, 
which implies 2h — f . In practice, we have calculated £ from 
EMD as a function of 2h for fractional Gaussian noises and 
motions with self-similarity exponent h, and have derived a 
relation ( = 0.94 (2h) + 0.1143. 

The variogram represent a second order structure func- 
tion: 

1 N ~* 

7fc = (JV-fc) H (s ' 1+fc ~ s " )2 ' (12) 

n— 1 

which scales with a time-lag k as 7^ = k 2h , h is denoted as 
selfsimilarity exponent, and s is a time series. Note that a 
Hurst exponent H > 1 implies that the process is a nonsta- 
tionary motion, and if the process is selfsimilar, the selfsim- 
ilarity exponent is h = H — 1. In our terminology a white 
noise process has Hurst exponent H = 0.5 and a Brownian 
motion has H = 1.5. 

From Figure 10a we observe three different scaling 
regimes for tSYM-H. For time scales less than a few hundred 
minutes it scales like an uncorrelated motion (h ~ 0.5). On 



time scales from a few hours to a week it scales as an an- 
tipersistent motion (h w 0.25 — 0.35 depending on analysis 
method), and on longer time scales than a week it is close 
to a stationary pink noise (h « 0). Similar behavior was ob- 
served for log AE in Rypdal and Rypdal [2010], but there the 
break between non-stationary motion and stationary noise 
(where h changes from h > to h w 0) occurs already after 
about 100 minutes, indicating the different time scales in- 
volved in ring current (storm) dynamics and electrojet (sub- 
storm) dynamics. 

Results for v indicate a regime with antipersistent motion 
(h = 0.25) up to a few hundred minutes, and then an uncor- 
related or weakly persistent motion (h = 0.5) up to a week. 
On longer time scales than this the variogram indicates that 
the process is stationary. 

The exponent h for B z can not be estimated from the 
variogram since it is difficult to obtain a linear fit to the con- 
cave curve in Figure 10b. The concavity is less pronounced 
in the curve derived from the EMD method in Figure 10a 
and C, = 0.47, corresponding to an antipersistent motion 
(h = 0.23), can be estimated on time scales up to a few 
hundred minutes. B z becomes stationary already after a 
few hundred minutes, which is similar to the behavior in 
log AE, as pointed out by Rypdal and Rypdal [2010]. In Ryp- 
dal and Rypdal [2011] the concavity of the variogram follows 
from modelling B z as a (multifractal) Ornstein-Uhlenbeck 
process with a strong damping term that confines the mo- 
tion on time scales longer than 100 minutes. Accounting for 




k (rain) 



Figure 10. a) The empirical variance estimates E ver- 
sus mean period T for each IF component in tSYM-H, 
B z , and v shown as log-log plots, b) The variogram 7*, 
shown in log-log plot. In both panels stars are for tSYM- 
H, diamonds for IMF B z , and triangles for v. Note that 
a generalization of the result £ = 2H — 2 in Flandrin and 
Gongalves [2004] yields 2h = (. 
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this confinement the "true" selfsimilarity exponent of the 
stochastic term turns out to be h m 0.5. Thus the antiper- 
sistence derived from the EMD analysis may be a spurious 
effect from this confinement. The conclusion in Rypdal and 
Rypdal [2011] is that B z and log AE behave as uncorrelated 
motions up to the scales of a few hours and become sta- 
tionary on scales longer than this. Moreover, the stochastic 
term modelling the two signals share the same multifractal 
spectrum. In comparison, tSYM-H and v are non-stationary 
motions on scales up to a week before they reach the sta- 
tionary regime. 

3.2. Change of determinism during storm times 

In Figure 11 we show L n and A(r) for tSYM-H and its 
surrogate time-series with randomized phases of Fourier co- 
efficients. We observe that L n and A(r) for the surrogate 
time series does not deviate from those computed from the 
original tSYM-H, indicating that the dynamics of tSYM-H 
is not low-dimensional and nonlinear. The same results are 
obtained for IMF B z and flow speed v (not shown here). 

In the following analysis we test for determinism in tSYM- 
H, B z and v for ten intense storms. The reference point in 
our analysis is the storm's main phase, and then we analyze 
all the data spanning the time interval three days before and 
three days after the storm in tSYM-H, B z and v. We com- 
pute L n for n = 6 with a time resolution of 12 hours. The 
choice of n — 6 from the L„ -curve is a compromise between 
clear separation between low-dimensional and stochastic dy- 
namics and small error bars (which increase with increasing 
n). In order to improve statistics for Lr,, we compute de- 
terminism using data from all ten storms. This means that 



each Lq is computed over 12 hours interval over 10 storms, 
which gives 12-60- 10 = 7200 points. As a reference, we com- 
pute 1/6 for the fO-U process, whose coefficients are fitted 
by the least square regression to the SYM-H index during 
investigated storms. In all computations, we use embedding 
dimension m = 7, time-delay r = 20, and 6=1. In Figure 
12a we plot L 6 for tSYM-H B z , v, and fO-U, and in Figure 
12b the D st index averaged over all ten storms is plotted, 
since this index shows precisely when the storm takes place. 
We can observe that Lq is essentially the same for B z , v and 
fO-U, and stays approximately constant during the course 
of a storm. However, Lg for tSYM-H increases during storm 
time. In order to demonstrate that the change in Lq is sig- 
nificant, we plot in Figure 13 L n for tSYM-H, where the 
triangles are the mean of L n computed 3 days before and 
after the storm for ten different storms. These curves rep- 
resent non-storm conditions. The upper curve (squares) is 
the mean over all ten storms computed at the time of the 
D 3 t minimum, i.e. it represents the L n -curve around storm 
onset. 

In addition, we test determ inism for the quantity SYM- 
H*=0.77 SYM-H- 11.9 y/Pdyn [Kozyra and Liemohn, 2003], 
where Pd yn is the Solar wind's dynamic pressure. SYM- 
H* is a corrected index where the effect of the magne- 
topause current due to Pd yn is subtracted, and thus rep- 
resents the ring-current contribution to SYM-H. In order 
to obtain stationary increments we analyze a transformed 
index; tSYM-H*=log(ci+c 2 SYM-H*), where ci = 1.7694, 
and C2 ~ 0.0292. Since some data points are missing in 
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Figure 11. a) L n : square symbols are for tSYM-H, and 
triangles are for this signal after randomization of phases 
of Fourier coefficients, b) A(r): squares is for tSYM-H, 
and triangles are after randomization of phases. 
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Figure 12. a) for IMF B z (squares), v (triangles), 
tSYM-H (diamonds), fO-U (stars) computed before, dur- 
ing and after storm onset. The values for L§ are com- 
puted using 12 hour intervals and are averaged over ten 
different storms, b) The D 3 t index averaged over the ten 
storms. 



TATJANA ZIVKOVIC, KRISTOFFER RYPDAL: MAGNETOSPHERE DURING STORMS 



X- 9 



Pdyn, we have made a linear interpolation over the missing 
points. It seems that the interpolation decreases determin- 
ism in tSYM-H* , and for reference we compute determinism 
for the interpolated tSYM-H, where the interpolation is done 
over the same points as in tSYM-H*, even though tSYM-H 
does not have missing points. In Figure 14 we plot L§ for 
tSYM-H*, tSYM-H, and ^J\P dyn ). We see that the deter- 
minism in tSYM-H* is lower than that in tSYM-H, but it 
still increases during the storm. On the other hand yJ\Pd yn ) 
shows no change in determinism during storm events. 

3.3. Discussion of results on determinism 

The determinism (as measured by Le) of the storm in- 
dex tSYM-H and tSYM-H* has been shown to exhibit a 
pronounced increase at storm time. A rather trivial expla- 
nation of this enhancement would be that it is caused by the 
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Figure 13. L n for tSYM-H, where the triangles are 
the mean of L n computed 3 days before and after the 
storm for ten different storms. These curves represent 
non-storm conditions. The upper curve (squares) is the 
mean over all ten storms computed at the time of the D st 
minimum, i.e. it represents the L n -curve around storm 
onset. Many curves are terminated for n max < 12 be- 
cause there were no boxes with more than n ma x passages 
of the phase-space trajectory. 
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Figure 14. Lg averaged over ten storms for time series 
where missing points have been interpolated; tSYM-H 
(diamonds), tSYM-H* (stars), -^/(Pdyn) (triangles) 



"trend" incurred by the wedge-shaped drop and recovery of 
the storm indices associated with a magnetic storm. We test 
this hypothesis by superposing such a wedge-shaped pulse 
to an fO-U process and compute Lq. Next, we take tSYM-H 
for ten storms and for each set of data subtract the wedge- 
shaped pulse (computed by a moving-average smoothing). 
The residual signal represents the "detrended" fluctuations. 
The result is shown in Figure 15 and reveals that the trend 
in fO-U process has no discernible influence on the determin- 
ism during the storm while, on the other hand, we observe 
that the enhancement of Lg around storm time persists in 
the "detrended" fluctuations. This result suggests that the 
increasing determinism during storms is a result of an en- 
hanced low-dimensional component in the storm indices. As 
mentioned in section 2.3 for low-dimensional dynamics, non- 
linearity may be important for the measure of determinism. 
For a nonlinear, low-dimensional system the destruction of 
nonlinear coupling by randomizing phases of Fourier coef- 
ficients will in general reduce the determinism, while for 
a linear, stochastic process we will observe no such effect. 
But what role will nonlinearity play if it is introduced in the 
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Figure 15. Lq: triangles are derived from an fO-U pro- 
cess with a "storm trend" imposed, diamonds are derived 
from the "detrended" tSYM-H. 
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Figure 16. The drift term in the fO-U equation com- 
puted from tSYM-H. The smooth solid curve is a six- 
order polynomial fit. 
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deterministic terms of a stochastic equation? The determin- 
istic term in the fractional Langevin equation representing 
the fO-U is a linear damping term. However, the best rep- 
resentation of the damping/drift term in an fO-U model for 
tSYM-H is not linear. Following Rypdal and Rypdal [2011], 
if y =tSYM-H, the drift term is given as the conditional 
probability density Sy(t, St) given that y(t) = y: 



M(y,St) = E[Sy(t,St)\y(t) = y]. 



(13) 



In fO-U M(y, St) is a linear function of y, but a polyno- 
mial fit to drift term derived from tSYM-H data requires a 
sixth order polynomial, confirming the nonlinearity of the 
tSYM-H process. This is shown in Figure 16. Next, we 
test determinism for the nonlinear fO-U process, whose scal- 
ing exponent h is estimated from the variogram of tSYM- 
H, and where m — 8 and r = 10 is used. Figure 17a 
shows L n for numerical realizations of this process compared 
with the same analysis after randomization of the phases of 
the Fourier coefficients. The result reveals that the non- 
linear fO-U process is not more deterministic than its ran- 
domized version. Next, we form the composite time series 
x = xl + 1.85x n , where xl is the solution of the Lorenz 
system and x n is the nonlinear fO-U process, both signals 
with zero mean and unit variance. Again, embedding di- 
mension m = 8 is used. Now the L„ -curve is lowered when 
the phases are randomized, as shown in Figure 17b, which 
confirms our conjecture that determinism is a measure of 
low-dimensionality. 
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Figure 17. L n . a) diamonds are derived from numerical 
solutions of the nonlinear fO-U. Triangles are from these 
solutions after randomization of phases of Fourier coeffi- 
cients, b) diamonds are derived from numerical solutions 
of the nonlinear fO-U with a solution of the x component 
of the Lorenz system superposed. Triangles are from the 
latter signals after randomization of phases. 



3.4. Change of predictability during magnetic storms 

Even though we deal with a predominantly stochastic sys- 
tem, its correlation and the degree of predictability changes 
in time, and our hypothesis is that abrupt transitions in the 
dynamics take place during events like magnetic storms and 
substorms. We therefore employ recurrence plot quantifica- 
tion analysis as a tool for detection of these transitions. 

We compute the average inverse diagonal line length 
r = (l^ 1 ) as defined in equation (11), but the same results 
can be drawn from other quantities that can be derived from 
the recurrence plot [Marwan et al, 2007]. T can be used as a 
proxy for the positive Lyapunov exponent in a system with 
chaotic dynamics, and is sensitive to the transition from reg- 
ular to chaotic behavior, as can be shown heuristically for 
the case of the Lorenz system, where we use a = 10, 6 = 8/3, 
and c is varied from 20 to 40, such that transient behavior 
is obtained. For c = 24.74 a Hopf bifurcation occurs, which 
corresponds to the onset of chaotic flows. In Figure 18a we 
plot a bifurcation diagram for the x component of the Lorenz 
system as a function of the parameter c, while in Figure 18b 
we show r for the x component again as a function of the 
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Figure 18. a) x component of the Lorenz system as a 
function of the parameter c. b) Y = (l^ 1 ) as a function 
of the parameter c. 
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Figure 19. a) V for tSYM-H (stars), tSYM-H* (circles), 
B z (squares), v (upward triangles), detrended tSYM-H 
(downward triangles) averaged over ten storms. Error 
bars represent standard deviation based on data from 
these ten storms. Time origin is defined by the minimum 
of the average D st index for the ten storms. 
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parameter c. Similar results have been obtained from the 
longest diagonal length, when applied to the logistic map 
[Trulla et al, 1996]. 

In the following analysis, we use embedding dimension 
m = 1, because the results do not seem dependent on m 
and because, in the case of stochastic or high-dimensional 
dynamics, a topological embedding cannot be achieved for 
any reasonable embedding dimension. This fact demon- 
strates the robustness of the recurrence-plot analysis, which 
responds to changes in the dynamics of the system even if 
it is a stochastic or high-dimensional system for which no 
proper phase-space reconstruction is possible. 

Since reduction in F means increase of predictability it 
may also be a signature of higher persistence in a stochastic 
signal. This motivates plotting F and 2h (computed as a 
linear fit from the variogram over the time scales up to 12 
hours) for solar wind parameters and magnetic indices. Fig- 
ure 19 shows T for tSYM-H, B z , v, tSYM-H* and detrended 
tSYM-H averaged over 10 magnetic storms. Figure 20 shows 
the same for 2h, but detrended tSYM-H is not shown since 
its 2h changes insignificantly during the course of the storm. 
We observe that the increase in the predictability and per- 
sistence does not occur simultaneously for all observables. 
While B z , tSYM-H, detrended tSYM-H and tSYM-H* get 
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Figure 20. a) 2h for tSYM-H (stars), tSYM-H* (circles), 
B z (squares) , v (triangles) averaged over the same storms 
as in the previous figure. 




the most predictable during or after the main phase of the 
storm, solar wind's flow speed becomes the most predictable 
prior to the storm's main phase. From a hundred realiza- 
tions of the fO-U process generated numerically with the 
coefficients in the stochastic equation fitted to model the 
tSYM-H signal, we find F — 0.4 and 2h — 1, in good agree- 
ment with the results obtained from the tSYM-H time se- 
ries. The general relationship between V and h can also be 
explored through numerical realizations of fO-U processes. 
Figure 21 shows F computed for varying h as a mean value 
of 100 realizations of such a process for each h. For persis- 
tent motions (h > 0.5) there is a linear dependence between 
F and h, and a best fit yields 



T « 0.72 - 0.57/i. 



(14) 



Figure 21. F vs. h computed from numerical realizations 
of the fO-U process. 



This analysis shows the importance of F as a universal 
measure for predictability: in low-dimensional systems it 
is a proxy for the Lyapunov exponent, while for persistent 
stochastic motions it is a measure of persistence through 
equation (14). 

4. Conclusions 

The storm index SYM-H and the solar wind observables 
(flow velocity v and IMF B z ) show no clear signatures of 
low-dimensional dynamics during quiet periods. However, 
low-dimensionality increases in SYM-H and SYM-H* during 
storm times, indicating that self-organization of the mag- 
netosphere takes place during magnetic storms. This con- 
clusion is drawn from the study of ten intense, magnetic 
storms in the period from 2000-2003. Even though our 
analysis shows no discernible change in determinism dur- 
ing magnetic storms for solar wind parameters, there is an 
enhancement of the predictability of the solar wind observ- 
ables as well as the geomagnetic storm indices during major 
storms. We interpret this as an increase in the persistence 
of the stochastic components of the signals. The increased 
persistence in the solar wind flow v, prior to the storm's 
main phase could indicate that v is more important driver 
than B z during magnetic storms. This is consistent with 
a reexamination of the solar wind-magnetosphere coupling 
functions done by Newell et al. [2006], who found that the 
most optimal function is of the form v 2 B sin 4 (6 , /2) 2//3 , where 
6 — arctan(_B !/ /_B z ). Also, it has been shown in Pulkkinen 
et al. [2007] through numerical simulation, that increased v 
changes the magnetospheric response from a steadily con- 
verting state to highly variable in both space and time. 

It has been shown in Nose et al. [2001] that the plasma 
sheet is the dominant source for the ring current based on 
the similarity in composition of the inner plasma sheet and 
ring current regions. During the main phase of the storm, 
ions from the plasma sheet are flowing to the inner mag- 
netosphere on the open drift paths and then move to the 
dayside magnetopause. In this storm phase the ring current 
is highly asymmetric, as was experimentally shown by ener- 
getic neutral atom imaging (see Kozyra and Liemohn [2003] 
and references therein). During the recovery phase, ions 
from the plasma sheet are trapped on closed drift paths, 
and form the symmetric ring current. Therefore, the in- 
crease in determinism of the ring-current (SYM-H*) during 
storms implies increased determinism in the plasma sheet as 
well. 

A magnetic storm is a coherent global phenomenon in- 
vesting a vast region of the inner magnetosphere, and imply- 
ing large scale correlation. The counterpart of this increase 
of coherence is the reduction of the spontaneous incoherent 
short time scale fluctuations. Consequently, one should ex- 
pect a reduction of the free degrees of freedom which implies 
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an increase of determinism, i.e. the possible emergence of a 
low-dimensionality. 

Analysis of predictability shows significant differences be- 
tween B z on one hand, and v and SYM-H on the other. 
While the former is a non-stationary, slightly anti-persistent 
motion up to time scales of approximately 100 minutes, and 
a pink noise on longer time scales, the latter are slightly per- 
sistent motions on scales up to several days and noises on 
longer time scales. These differences indicate the different 
role the solar wind B z and the velocity v play in driving 
the substorm and storm current systems; B z is important 
in substorm dynamics which will be studied in a separate 
paper, while v is a major driver of storms. 
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